rm(list = ls())
require(zTree)
require(stargazer)
require(mfx)
require(stats4)
require(nlWaldTest)
require("multiwayvcov")
require(msm)
require(miceadds)
require(restriktor)

zTT1= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/August2Expansionpath2021.xls")
subjects1=zTT1$subjects

zTT2= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/September21ExpansionPath2021.xls")
subjects2=zTT2$subjects
subjects2$Subject=subjects2$Subject+100

zTT3= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/September 29ExpansionPath2021.xls")
subjects3=zTT3$subjects
subjects3$Subject=subjects3$Subject+200

zTT4= zTreeTables("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/Ambueletalexperiments/Data/October13expansionpath2021.xls")
subjects4=zTT4$subjects
subjects4$Subject=subjects4$Subject+300


subjects=rbind(subjects1,subjects2,subjects3,subjects4)

#Clean data
subjects$responder=(subjects$"response1[50]">0)
subjects=subset(subjects,responder>0)

uidlist=unique(subjects$Subject)
incentivec=c(5,95)

#create a frame with observations
obsframe=data.frame(uid=vector('numeric'),block=vector('numeric'), incentive=vector('numeric'), number=vector('numeric'), state=vector('numeric'), response=vector('numeric'),react=vector('numeric'))
for(i in 1:length(uidlist)){
playframe=subjects[i,]
uidplace=uidlist[i]

#block 1
for (j in 1:50){
blockplace=1
incentiveplace=incentivec[playframe[1,683]]
numberplace=j
stateplace=playframe[1,(736+j)]
responseplace=playframe[1,(72+j)]
reactplace=playframe[1,(122+j)]
addframe=data.frame(uid=uidplace,block=blockplace, incentive=incentiveplace, number=numberplace, state=stateplace, response=responseplace,react=reactplace)
obsframe=rbind(obsframe,addframe)
}

#block 2
for (j in 1:50){i
blockplace=2
incentiveplace=incentivec[playframe[1,684]]
numberplace=j
stateplace=playframe[1,(786+j)]
responseplace=playframe[1,(222+j)]
reactplace=playframe[1,(272+j)]
addframe=data.frame(uid=uidplace,block=blockplace, incentive=incentiveplace, number=numberplace, state=stateplace, response=responseplace,react=reactplace)
obsframe=rbind(obsframe,addframe)

}

}



#define correctness
obsframe$correct=(obsframe$state==obsframe$response)

#mean correctness
correct5=sum(obsframe$correct*(obsframe$incentive==5))/sum(obsframe$incentive==5)
correct95=sum(obsframe$correct*(obsframe$incentive==95))/sum(obsframe$incentive==95)
correctness=c(correct5,correct95)

obsframe5=subset(obsframe,incentive==5)
obsframe95=subset(obsframe,incentive==95)

#change terminology to work with older scripts
expframe=obsframe

expframe$incentive5=(expframe$incentive==5)
expframe$incentive95=(expframe$incentive==95)

expframe$Bresponse <-expframe$response==2
expframe$Bstate <- expframe$state==2
expframe$Aresponse <-expframe$response==1
expframe$Astate <- expframe$state==1

expframe5=subset(expframe,incentive==5)
expframe95=subset(expframe,incentive==95)

#NIAS in aggregate

#NIAS in aggregate
#Base Regression
expframe$incentivedec=expframe$incentive/100
modelNIAS=lm(Aresponse~factor(Astate+incentivedec)+0,data=expframe)

#Adjust COV
vcovmodNIAScluster=cluster.vcov(modelNIAS, expframe$uid)

modelniassterr=vector("numeric")
#Model Nias sterrors
for(i in 1:4){
modelniassterr[i]=sqrt(vcovmodNIAScluster[i,i])
}


#Check NIAS and NIAC on the individual level
#NIAS


NIASAframe=data.frame(uid=vector(),A_state_avg_effect=vector(),Pvalue=vector(),incentive=vector())
for (i in 1:length(uidlist)){
for(j in 1:2){
addframe=0
playframe=subset(expframe,uid==uidlist[i] & incentive==incentivec[j])
if(sum(playframe$Aresponse)==0|sum(playframe$Bresponse)==0){
addframe=data.frame(uid=uidlist[i],A_state_coefficient=0,Pvalue=1,incentive=incentivec[j])
NIASAframe=rbind(NIASAframe,addframe)
}else{
model=logitmfx(Aresponse~Astate,data=playframe,robust=TRUE) 
addframe=data.frame(uid=uidlist[i],A_state_coefficient=model[[1]][1],Pvalue=model[[1]][4],incentive=incentivec[j])
NIASAframe=rbind(NIASAframe,addframe)
}
}
}

#Check NIAC 

restrictmat=matrix(c(0,1),nrow=1,ncol=2, byrow=TRUE)

NIACindiframe=data.frame(uid=vector('numeric'),fstat=vector('numeric'),pvalue=vector('numeric'))
for(k in 1:length(uidlist)){
playframe=subset(expframe,uid==uidlist[k])

playmodel=glm(correct~incentive95,data=playframe,family=binomial())
playtest=conTest(playmodel,constraints=restrictmat,type="B",rhs=c(0))
addframe=data.frame(uid=uidlist[k],fstat=playtest$Ts,pvalue=playtest$pvalue)

NIACindiframe=rbind(NIACindiframe,addframe)
}
sum(NIACindiframe$pvalue<0.05)/length(NIACindiframe$pvalue)

#NIAS and NIAC together
comboframe=data.frame(uid=vector('numeric'),NIASviol=vector('numeric'),NIACviol=vector('numeric'))
for(i in 1:length(uidlist)){
NIACplayframe=subset(NIACindiframe,uid==uidlist[i])
NIASplayframe=subset( NIASAframe,uid==uidlist[i])
addframe=data.frame(uid=uidlist[i],NIASviol=(sum((NIASplayframe$A_state_coefficient<0)*(NIASplayframe$Pvalue<0.05))>=1),NIACviol=(NIACplayframe$pvalue<0.05))
comboframe=rbind(comboframe,addframe)
}

comboframeclean=comboframe
NIASonly=100*sum((comboframeclean$NIASviol==1)*(comboframeclean$NIACviol==FALSE))/length(uidlist)
NIAConly=100*sum((comboframeclean$NIASviol==0)*(comboframeclean$NIACviol==TRUE))/length(uidlist)
both=100*sum((comboframeclean$NIASviol==1)*(comboframeclean$NIACviol==TRUE))/length(uidlist)
neither=100*sum((comboframeclean$NIASviol==0)*(comboframeclean$NIACviol==FALSE))/length(uidlist)
violatereport=data.frame(violate=c("NIAS only","NIAC only","Both","Neither"),percent=c(NIASonly,NIAConly,both,neither))

#Table A5.3************************************
stargazer(violatereport,summary=FALSE,rownames=FALSE)
#************************************


